Revised Interoceptive Accuracy Scale (IAS-R)

Analysis

Code

library(tidyverse)
library(easystats)
library(patchwork)
library(lavaan)
library(ggraph)
library(tidySEM)
library(EGAnet)

Study 1 - Reanalysis of Murphy et al. (2020)

Reanalyze their data to confirm the factor structure of the IAS.

Data Preparation

Code

# https://osf.io/3m5nh
df1 <- haven::read_sav("../data/Murphy2020/Study 1.sav") |> 
  mutate_all(as.numeric) |> 
  mutate(Gender = as.character(ifelse(Gender == 1, "Male", ifelse(Gender == 2, "Female", "Other")))) 

# https://osf.io/3m5nh
df2 <- haven::read_sav("../data/Murphy2020/Study 6 IAS.sav") |> 
  mutate_all(as.numeric) |> 
  mutate(Gender = as.character(ifelse(Gender == 1, "Male", ifelse(Gender == 2, "Female", "Other")))) 

# Campos (2022) - Study 1 (https://osf.io/j6ef3)
df3 <- haven::read_sav("../data/Campos2022/Dataset_Test.sav") |> 
  mutate_all(as.numeric) |> 
  mutate(Gender = as.character(ifelse(Sex == 1, "Male", ifelse(Sex == 0, "Female", "Other")))) |> 
  rename(Heart=IAS1, Hungry=IAS2, Breathing=IAS3, Thirsty=IAS4,
         Urinate=IAS5, Defecate=IAS6, Taste=IAS7, Vomit=IAS8, Sneeze=IAS9,
         Cough=IAS10, Temp=IAS11, Sex_arousal=IAS12, Wind=IAS13, Burp=IAS14,
         Muscles=IAS15, Bruise=IAS16, Pain=IAS17, Blood_Sugar=IAS18,
         Affective_touch=IAS19, Tickle=IAS20, Itch=IAS21)

# https://osf.io/3eztd
df4 <- haven::read_sav("../data/Lin2023/Study 1 & 3.sav") |> 
  mutate(Gender = as.character(ifelse(sex_dummy == 1, "Male", ifelse(sex_dummy == 0, "Female", "Other")))) |> 
  rename(Age=age, 
         Heart=Heart, Hungry=HUNGR, Breathing=BREAT, Thirsty=Thirs,
         Urinate=URINA, Defecate=Defec, Taste=TASTE, Vomit=VOMIT, Sneeze=Sneez,
         Cough=COUGH, Temp=TEMPE, Sex_arousal=SEXAR, Wind=WIND, Burp=Burp,
         Muscles=MUSCL, Bruise=Bruis, Pain=PAIN, Blood_Sugar=BloSu,
         Affective_touch=Touch, Itch=ITCH)  # No tickle because same Chinese character

# https://osf.io/3eztd
df5 <- haven::read_sav("../data/Lin2023/Study 2.sav") |> 
  mutate(Gender = as.character(ifelse(Sex == "男", "Male", ifelse(Sex == "女", "Female", "Other")))) |> 
  rename(Heart=Heart, Hungry=HUNGR, Breathing=BREAT, Thirsty=Thirs,
         Urinate=URINA, Defecate=Defec, Taste=TASTE, Vomit=VOMIT, Sneeze=Sneez,
         Cough=COUGH, Temp=TEMPE, Sex_arousal=SEXAR, Wind=WIND, Burp=Burp,
         Muscles=MUSCL, Bruise=Bruis, Pain=PAIN, Blood_Sugar=BloSu,
         Affective_touch=Touch, Itch=ITCH)  # No tickle because same Chinese character

df6 <- read.csv("https://raw.githubusercontent.com/DominiqueMakowski/PHQ4R/main/study2/data/data.csv") |>
  rename(Heart=IAS_1, Hungry=IAS_2, Breathing=IAS_3, Thirsty=IAS_4, 
         Urinate=IAS_5, Defecate=IAS_6, Taste=IAS_7, Vomit=IAS_8, Sneeze=IAS_9, 
         Cough=IAS_10, Temp=IAS_11, Sex_arousal=IAS_12, Wind=IAS_13, Burp=IAS_14,
         Muscles=IAS_15, Bruise=IAS_16, Pain=IAS_17, Blood_Sugar=IAS_18,
         Affective_touch=IAS_19, Tickle=IAS_20, Itch=IAS_21)

df7 <- read.csv("https://raw.githubusercontent.com/RealityBending/IllusionGameReliability/main/data/preprocessed_questionnaires.csv") |>
  rename(Gender=Sex, Heart=Item_IAS_Interoception_1, Hungry=Item_IAS_Interoception_2, 
         Breathing=Item_IAS_Interoception_3, Thirsty=Item_IAS_Interoception_4, 
         Temp=Item_IAS_Interoception_5, Sex_arousal=Item_IAS_Interoception_6,
         Urinate=Item_IAS_Elimination_1, Defecate=Item_IAS_Elimination_2,
         Vomit=Item_IAS_Elimination_3, Wind=Item_IAS_Expulsion_1, 
         Burp=Item_IAS_Expulsion_2, Sneeze=Item_IAS_Expulsion_3,
         Muscles=Item_IAS_Nociception_1, Bruise=Item_IAS_Nociception_2,
         Pain=Item_IAS_Nociception_3, Affective_touch=Item_IAS_Skin_1, 
         Tickle=Item_IAS_Skin_2, Itch=Item_IAS_Skin_3) |>
  filter(!is.na(Urinate))

Participants

• Sample 1: Data from Murphy’s (2020) study 1, downloaded from OSF, included 451 participants (Mean age = 25.8, SD = 8.4, range: [18, 69]; Gender: 69.4% women, 29.5% men, 1.11% non-binary).
• Sample 2: Data from Murphy’s (2020) study 6, downloaded from OSF, included 375 participants (Mean age = 35.3, SD = 16.9, range: [18, 91]; Gender: 70.1% women, 28.5% men, 1.33% non-binary).
• Sample 3: Data from Campos’ (2022) study 1, downloaded from OSF, included 515 participants (Mean age = 30.7, SD = 10.5, range: [18, 72]; Sex: 0.0% females, 0.0% males, 100.0% other; Gender: 59.6% women, 40.4% men, 0.00% non-binary).
• Sample 4: Data from Lin’s (2023) study 1 and 3, downloaded from OSF, included 1166 participants (Mean age = 32.5, SD = 8.4, range: [16, 60]; Gender: 57.0% women, 43.0% men, 0.00% non-binary).
• Sample 5: Data from Lin’s (2023) study 2, downloaded from OSF, included 500 participants (Mean age = 37.4, SD = 7.4, range: [20, 60]; Sex: 0.0% females, 0.0% males, 100.0% other; Gender: 56.2% women, 43.8% men, 0.00% non-binary).
• Sample 6: New dataset from Makowski (under preparation), included 251 participants (Mean age = 37.5, SD = 13.2, range: [17, 76]; Gender: 53.0% women, 42.6% men, 2.79% non-binary, 1.59% missing; Education: Bachelor, 39.84%; Doctorate, 4.78%; High School, 31.47%; Master, 19.12%; missing, 0.40%; Other, 4.38%).
• Sample 7: New dataset from Makowski (under preparation), included 485 participants (Mean age = 30.1, SD = 10.1, range: [18, 73]; Gender: 50.3% women, 49.7% men, 0.00% non-binary; Education: Bachelor, 45.15%; Doctorate, 1.86%; High school, 34.43%; Master, 15.88%; Other, 2.47%; Prefer not to say, 0.21%).

Total N = 3743.

Descriptive

Distribution

Code

vars <- c("Heart", "Hungry", "Breathing", "Thirsty", "Urinate", "Defecate", "Taste", "Vomit", "Sneeze", "Cough", "Temp",
          "Sex_arousal", "Wind", "Burp", "Muscles", "Bruise", "Pain", "Blood_Sugar", "Affective_touch", "Tickle", "Itch")
  
dens1 <- estimate_density(select(df1, all_of(vars)), method="kernSmooth") |> 
  mutate(Sample = "Sample 1")
dens2 <- estimate_density(select(df2, all_of(vars)), method="kernSmooth") |> 
  mutate(Sample = "Sample 2")
dens3 <- estimate_density(select(df3, all_of(vars)), method="kernSmooth") |> 
  mutate(Sample = "Sample 3")
dens4 <- estimate_density(select(df4, all_of(setdiff(vars, "Tickle"))), method="kernSmooth") |> 
  mutate(Sample = "Sample 4")
dens5 <- estimate_density(select(df5, all_of(setdiff(vars, "Tickle"))), method="kernSmooth") |> 
  mutate(Sample = "Sample 5")
dens6 <- estimate_density(select(df6, all_of(vars)), method="kernSmooth") |> 
  mutate(Sample = "Sample 6")
dens7 <- estimate_density(select(df7, all_of(setdiff(vars, c("Taste", "Cough", "Blood_Sugar")))), method="kernSmooth") |> 
  mutate(Sample = "Sample 7")

rbind(dens1, dens2, dens3, dens4, dens5, dens6, dens7) |> 
  ggplot(aes(x = x, y = y, color = Parameter)) +
  geom_line() +
  theme_minimal() +
  labs(title = "Distribution of Responses", x = "Response", y = "Density", color = "Item") +
  facet_wrap(~Sample, scales = "free", nrow=4)

Code

data1 <- normalize(select(df1, all_of(dens1$Parameter)), range=c(1, 5))
data2 <- normalize(select(df2, all_of(dens2$Parameter)), range=c(1, 5))
data3 <- normalize(select(df3, all_of(dens3$Parameter)), range=c(1, 5))
data4 <- normalize(select(df4, all_of(dens4$Parameter)), range=c(1, 5))
data5 <- normalize(select(df5, all_of(dens5$Parameter)), range=c(1, 5))
data6 <- select(df6, all_of(dens6$Parameter))
data7 <- select(df7, all_of(dens7$Parameter))

data_all <- rbind(data1, data2, data3, 
                  mutate(data4, Tickle=NA), mutate(data5, Tickle=NA), 
                  data6, mutate(data7, Taste=NA, Cough=NA, Blood_Sugar=NA)) 

Correlations

An overall positive intercorrelation pattern, with no clear structure emerging.

All data

Code

make_correlation <- function(df) {
  correlation::correlation(df, redundant=TRUE) |> 
    correlation::cor_sort() |> 
    correlation::cor_lower() |> 
    mutate(val = paste0(insight::format_value(r), format_p(p, stars_only=TRUE))) |>
    mutate(Parameter2 = fct_rev(Parameter2)) |>
    ggplot(aes(x=Parameter1, y=Parameter2)) +
    geom_tile(aes(fill = r), color = "white") +
    geom_text(aes(label = val), size = 3) +
    labs(title = "Correlation Matrix") +
    scale_fill_metro_c(limit = c(0, 0.75), guide = guide_colourbar(ticks=FALSE)) +
    theme_minimal() +
    theme(legend.title = element_blank(),
          axis.title.x = element_blank(),
          axis.title.y = element_blank(),
          axis.text.x = element_text(angle = 45, hjust = 1))
}

make_correlation(data_all)

Sample 1

Code

make_correlation(data1)

Sample 2

Code

make_correlation(data2)

Sample 3

Code

make_correlation(data3)

Sample 4

Code

make_correlation(data4)

Sample 5

Code

make_correlation(data5)

Sample 6

Code

make_correlation(data6)

Sample 7

Code

make_correlation(data7)

EGA

See https://r-ega.net/articles/ega.html for details.

Unique Variable Analysis (UVA)

All data

Code

uva <- EGAnet::UVA(data = data_all, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.306

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i   node_j   wto
  Hungry  Thirsty 0.222
 Urinate Defecate 0.219
  Sneeze    Cough 0.204

Code

uva$keep_remove

$keep
[1] "Itch"

$remove
[1] "Tickle"

Sample 1

Code

uva <- EGAnet::UVA(data = data1, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.286
   Wind   Burp 0.275

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i   node_j   wto
  Sneeze    Cough 0.244
 Urinate Defecate 0.241

Code

uva$keep_remove

NULL

Sample 2

Code

uva <- EGAnet::UVA(data = data2, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.350
 Sneeze  Cough 0.317
   Wind   Burp 0.309

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i   node_j   wto
 Urinate Defecate 0.278

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

Code

uva$keep_remove

$keep
[1] "Cough" "Burp"  "Itch" 

$remove
[1] "Sneeze" "Wind"   "Tickle"

Sample 3

Code

uva <- EGAnet::UVA(data = data3, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.445
 Sneeze  Cough 0.318

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i node_j   wto
   Wind   Burp 0.256

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i      node_j   wto
  Bruise Blood_Sugar 0.219
 Urinate    Defecate 0.217
   Heart   Breathing 0.208

Code

uva$keep_remove

$keep
[1] "Sneeze" "Itch"  

$remove
[1] "Cough"  "Tickle"

Sample 4

Code

uva <- EGAnet::UVA(data = data4, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i    node_j   wto
 Urinate  Defecate 0.237
   Heart Breathing 0.218
  Hungry   Thirsty 0.213

Code

uva$keep_remove

NULL

Sample 5

Code

uva <- EGAnet::UVA(data = data5, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i    node_j   wto
 Urinate  Defecate 0.277
   Heart Breathing 0.267

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i node_j   wto
   Wind   Burp 0.206

Code

uva$keep_remove

NULL

Sample 6

Code

uva <- EGAnet::UVA(data = data6, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

Code

uva$keep_remove

NULL

Sample 7

Code

uva <- EGAnet::UVA(data = data7, cut.off = 0.3)
uva

Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.297

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i    node_j   wto
  Heart Breathing 0.209

Code

uva$keep_remove

NULL

Item Stability

Code

plots <- list()
for(model in c("glasso", "TMFG")) {
  for(algo in c("walktrap", "louvain")) {
    for(type in c("ega", "ega.fit")) {  # "hierega"
      print(type)
      ega <- EGAnet::bootEGA(
        data = data_all,
        seed=123,
        model=model,
        algorithm=algo,
        EGA.type=type,
        type="resampling",
        plot.itemStability=FALSE,
        verbose=FALSE)
      plots[[length(plots) + 1]] <- plot(ega) + 
        labs(title = paste0(model, " (", algo, ") - ", type))
      }
    }
}

[1] "ega"
[1] "ega.fit"
[1] "ega"
[1] "ega.fit"
[1] "ega"
[1] "ega.fit"
[1] "ega"
[1] "ega.fit"

Code

length(plots)

[1] 8

Code

patchwork::wrap_plots(plots, nrow = 4)

Exploratory Factor Analysis (EFA)

How Many Factors

All data

Code

n <- parameters::n_factors(data_all, n_max = 12)

plot(n) 

Code

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 4, 5 ), ], row.names = FALSE)

n_Factors	Method	Family
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
1	Velicer’s MAP	Velicers_MAP
4	beta	Multiple_regression
4	Scree (SE)	Scree_SE
5	Optimal coordinates	Scree
5	Parallel analysis	Scree
5	Kaiser criterion	Scree

Sample 1

Code

n <- parameters::n_factors(data1, n_max = 12)

plot(n)

Code

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 4), ], row.names = FALSE)

n_Factors	Method	Family
1	Bentler	Bentler
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
1	Velicer’s MAP	Velicers_MAP
4	beta	Multiple_regression
4	Optimal coordinates	Scree
4	Parallel analysis	Scree
4	Kaiser criterion	Scree

Sample 2

Code

n <- parameters::n_factors(data2, n_max = 12)

plot(n)

Code

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 4), ], row.names = FALSE)

n_Factors	Method	Family
1	Bentler	Bentler
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
4	beta	Multiple_regression
4	Optimal coordinates	Scree
4	Parallel analysis	Scree
4	Kaiser criterion	Scree

Sample 3

Code

n <- parameters::n_factors(data3, n_max = 12)

plot(n) 

Code

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 4), ], row.names = FALSE)

n_Factors	Method	Family
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
4	beta	Multiple_regression
4	Optimal coordinates	Scree
4	Parallel analysis	Scree
4	Kaiser criterion	Scree
4	Scree (SE)	Scree_SE
4	BIC	BIC

Sample 4

Code

n <- parameters::n_factors(data4, n_max = 12)

plot(n) 

Code

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 4), ], row.names = FALSE)

n_Factors	Method	Family
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
1	Velicer’s MAP	Velicers_MAP
4	beta	Multiple_regression
4	Scree (SE)	Scree_SE

Sample 5

Code

n <- parameters::n_factors(data5, n_max = 12)

plot(n) 

Code

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 4), ], row.names = FALSE)

n_Factors	Method	Family
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
1	Velicer’s MAP	Velicers_MAP
4	beta	Multiple_regression
4	BIC	BIC

Sample 6

Code

n <- parameters::n_factors(data6, n_max = 12)

plot(n) 

Code

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 4), ], row.names = FALSE)

n_Factors	Method	Family
1	Bentler	Bentler
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
1	Velicer’s MAP	Velicers_MAP
4	beta	Multiple_regression
4	Optimal coordinates	Scree
4	Parallel analysis	Scree
4	Kaiser criterion	Scree
4	Scree (SE)	Scree_SE

Sample 7

Code

n <- parameters::n_factors(data7, n_max = 12)

plot(n) 

Code

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 4), ], row.names = FALSE)

n_Factors	Method	Family
1	Acceleration factor	Scree
1	Scree (R2)	Scree_SE
1	VSS complexity 1	VSS
1	Velicer’s MAP	Velicers_MAP
4	beta	Multiple_regression
4	Optimal coordinates	Scree
4	Parallel analysis	Scree
4	Kaiser criterion	Scree
4	BIC	BIC
4	BIC (adjusted)	BIC

EFA Models

All data

Code

efa5 <- parameters::factor_analysis(data_all, n=5, rotation = "oblimin", sort = TRUE)

Loading required namespace: GPArotation

Code

plot(efa5)

Code

display(efa5)

Rotated loadings from Factor Analysis (oblimin-rotation)

Variable	MR3	MR2	MR1	MR5	MR4	Complexity	Uniqueness
Cough	0.72	-5.89e-04	-0.06	-5.43e-03	0.03	1.02	0.50
Sneeze	0.67	-0.03	0.04	0.06	-0.05	1.04	0.53
Burp	0.55	0.13	0.06	-0.03	-0.02	1.15	0.59
Wind	0.48	0.11	-9.06e-03	-1.90e-03	0.06	1.15	0.68
Vomit	0.35	0.03	0.24	0.14	-0.09	2.32	0.68
Itch	0.05	0.71	-0.10	-0.06	0.12	1.12	0.46
Tickle	0.08	0.68	-0.03	0.06	4.30e-03	1.05	0.45
Bruise	8.13e-04	0.45	0.22	0.06	-0.08	1.56	0.68
Blood_Sugar	0.06	0.36	0.14	0.04	-0.06	1.43	0.79
Affective_touch	-0.05	0.36	0.34	0.18	-0.16	2.91	0.64
Muscles	0.03	0.25	0.21	0.17	0.07	2.95	0.72
Taste	0.16	0.18	0.17	0.13	0.05	3.96	0.75
Urinate	0.10	-0.08	0.55	0.01	0.20	1.39	0.55
Defecate	0.22	-0.04	0.46	-0.01	0.15	1.69	0.60
Pain	-0.01	0.22	0.28	0.11	0.23	3.22	0.65
Heart	1.28e-03	-0.03	-0.01	0.66	-0.01	1.01	0.59
Breathing	0.07	7.82e-03	-0.08	0.56	0.14	1.20	0.61
Sex_arousal	0.10	0.07	0.18	0.25	0.12	2.87	0.73
Thirsty	-1.95e-04	0.05	0.11	6.26e-03	0.62	1.07	0.53
Hungry	0.01	0.04	-5.98e-03	0.12	0.60	1.08	0.55
Temp	0.11	0.12	0.14	0.15	0.27	3.12	0.68

The 5 latent factors (oblimin rotation) accounted for 38.36% of the total variance of the original data (MR3 = 10.15%, MR2 = 9.46%, MR1 = 6.49%, MR5 = 6.20%, MR4 = 6.05%).

Sample 1

Code

efa4 <- parameters::factor_analysis(data2, n=4, rotation = "oblimin", sort = TRUE)
plot(efa4)

Code

display(efa4)

Rotated loadings from Factor Analysis (oblimin-rotation)

Variable	MR1	MR2	MR3	MR4	Complexity	Uniqueness
Urinate	0.69	-0.18	0.13	-0.04	1.22	0.54
Defecate	0.68	-0.13	0.04	0.06	1.10	0.55
Breathing	0.56	0.03	0.06	7.96e-03	1.03	0.63
Hungry	0.55	0.15	-7.15e-03	0.02	1.16	0.60
Thirsty	0.52	0.05	0.14	-0.01	1.16	0.62
Sex_arousal	0.50	0.14	-0.04	0.20	1.49	0.56
Temp	0.46	0.23	0.16	-0.01	1.77	0.52
Pain	0.46	0.43	-0.05	-0.03	2.03	0.49
Taste	0.40	0.21	-0.07	0.12	1.77	0.69
Heart	0.40	0.06	-0.04	0.15	1.36	0.76
Muscles	0.35	0.24	0.07	0.09	2.07	0.65
Vomit	0.32	0.06	0.23	0.11	2.20	0.68
Tickle	0.02	0.69	0.01	0.04	1.01	0.47
Itch	-0.02	0.67	0.18	0.03	1.16	0.41
Bruise	0.07	0.44	0.15	0.06	1.32	0.65
Blood_Sugar	-0.08	0.38	0.16	0.09	1.61	0.77
Affective_touch	0.38	0.38	-0.13	0.05	2.26	0.63
Sneeze	0.05	0.05	0.71	0.01	1.02	0.41
Cough	0.04	0.02	0.71	0.04	1.01	0.43
Wind	7.52e-03	-0.04	-5.18e-03	0.94	1.00	0.15
Burp	-0.01	0.27	0.17	0.47	1.91	0.48

The 4 latent factors (oblimin rotation) accounted for 44.25% of the total variance of the original data (MR1 = 17.71%, MR2 = 11.72%, MR3 = 7.65%, MR4 = 7.17%).

Sample 2

Code

efa4 <- parameters::factor_analysis(data2, n=4, rotation = "oblimin", sort = TRUE)
plot(efa4)

Code

display(efa4)

Rotated loadings from Factor Analysis (oblimin-rotation)

Variable	MR1	MR2	MR3	MR4	Complexity	Uniqueness
Urinate	0.69	-0.18	0.13	-0.04	1.22	0.54
Defecate	0.68	-0.13	0.04	0.06	1.10	0.55
Breathing	0.56	0.03	0.06	7.96e-03	1.03	0.63
Hungry	0.55	0.15	-7.15e-03	0.02	1.16	0.60
Thirsty	0.52	0.05	0.14	-0.01	1.16	0.62
Sex_arousal	0.50	0.14	-0.04	0.20	1.49	0.56
Temp	0.46	0.23	0.16	-0.01	1.77	0.52
Pain	0.46	0.43	-0.05	-0.03	2.03	0.49
Taste	0.40	0.21	-0.07	0.12	1.77	0.69
Heart	0.40	0.06	-0.04	0.15	1.36	0.76
Muscles	0.35	0.24	0.07	0.09	2.07	0.65
Vomit	0.32	0.06	0.23	0.11	2.20	0.68
Tickle	0.02	0.69	0.01	0.04	1.01	0.47
Itch	-0.02	0.67	0.18	0.03	1.16	0.41
Bruise	0.07	0.44	0.15	0.06	1.32	0.65
Blood_Sugar	-0.08	0.38	0.16	0.09	1.61	0.77
Affective_touch	0.38	0.38	-0.13	0.05	2.26	0.63
Sneeze	0.05	0.05	0.71	0.01	1.02	0.41
Cough	0.04	0.02	0.71	0.04	1.01	0.43
Wind	7.52e-03	-0.04	-5.18e-03	0.94	1.00	0.15
Burp	-0.01	0.27	0.17	0.47	1.91	0.48

The 4 latent factors (oblimin rotation) accounted for 44.25% of the total variance of the original data (MR1 = 17.71%, MR2 = 11.72%, MR3 = 7.65%, MR4 = 7.17%).

Sample 3

Code

efa4 <- parameters::factor_analysis(data3, n=4, rotation = "oblimin", sort = TRUE)
plot(efa4)

Code

display(efa4)

Rotated loadings from Factor Analysis (oblimin-rotation)

Variable	MR1	MR3	MR2	MR4	Complexity	Uniqueness
Cough	0.88	-0.08	-0.04	0.11	1.05	0.29
Sneeze	0.74	0.05	1.42e-03	-0.03	1.01	0.41
Vomit	0.55	0.12	0.04	-0.09	1.17	0.58
Burp	0.55	0.08	0.22	-0.13	1.47	0.48
Wind	0.53	0.10	0.21	-0.25	1.84	0.48
Temp	0.46	0.34	-0.06	0.09	1.97	0.52
Taste	0.29	0.18	0.10	0.12	2.37	0.74
Hungry	-8.85e-03	0.70	-0.01	-0.09	1.03	0.54
Urinate	0.08	0.66	-0.12	0.11	1.15	0.52
Thirsty	-0.08	0.64	-0.03	0.14	1.14	0.62
Defecate	0.10	0.63	-3.50e-05	-0.07	1.08	0.53
Breathing	-0.10	0.49	0.28	-0.07	1.76	0.66
Sex_arousal	0.17	0.44	0.17	-0.22	2.20	0.58
Pain	0.21	0.41	0.08	0.18	2.00	0.58
Heart	0.05	0.28	0.24	0.08	2.23	0.77
Muscles	0.24	0.27	0.20	0.19	3.66	0.61
Tickle	0.02	-0.05	0.84	-0.04	1.01	0.31
Itch	0.03	0.02	0.74	0.11	1.05	0.39
Affective_touch	0.08	0.23	0.36	0.12	2.08	0.67
Bruise	0.17	0.07	0.29	0.45	2.08	0.55
Blood_Sugar	4.11e-03	0.11	0.35	0.43	2.07	0.60

The 4 latent factors (oblimin rotation) accounted for 45.69% of the total variance of the original data (MR1 = 15.50%, MR3 = 15.46%, MR2 = 11.02%, MR4 = 3.71%).

Sample 4

Code

efa5 <- parameters::factor_analysis(data4, n=5, rotation = "oblimin", sort = TRUE)
plot(efa5)

Code

display(efa5)

Rotated loadings from Factor Analysis (oblimin-rotation)

Variable	MR1	MR3	MR5	MR4	MR2	Complexity	Uniqueness
Cough	0.67	-0.07	8.90e-03	0.06	0.03	1.04	0.54
Sneeze	0.63	-0.04	0.11	0.05	-0.08	1.11	0.54
Burp	0.58	0.18	-0.03	-0.09	0.03	1.25	0.59
Wind	0.54	0.08	-0.06	-0.03	0.11	1.17	0.67
Vomit	0.33	0.06	0.21	0.15	-0.19	3.02	0.70
Taste	0.24	0.19	-4.60e-03	0.06	0.18	3.00	0.77
Bruise	-0.03	0.63	0.05	0.05	-0.03	1.03	0.58
Affective_touch	0.07	0.53	0.03	-0.01	-0.04	1.05	0.67
Blood_Sugar	0.07	0.43	0.15	0.10	-0.18	1.82	0.67
Itch	0.11	0.41	-0.01	0.05	0.15	1.45	0.71
Muscles	0.07	0.40	-0.05	0.02	0.19	1.56	0.75
Urinate	-0.04	0.06	0.64	6.75e-03	0.06	1.04	0.55
Defecate	0.11	0.02	0.59	-0.03	0.02	1.08	0.57
Heart	-0.04	0.04	-0.03	0.66	-0.04	1.02	0.59
Breathing	0.09	-2.40e-03	-0.02	0.51	0.12	1.17	0.65
Sex_arousal	0.07	0.18	0.07	0.24	0.13	2.86	0.78
Thirsty	0.05	-0.05	0.24	0.03	0.49	1.51	0.60
Hungry	0.09	-0.10	0.15	0.21	0.38	2.25	0.67
Temp	0.09	0.17	7.51e-03	0.11	0.36	1.81	0.72
Pain	-3.56e-03	0.30	0.14	0.04	0.33	2.38	0.67

The 5 latent factors (oblimin rotation) accounted for 35.08% of the total variance of the original data (MR1 = 10.08%, MR3 = 8.61%, MR5 = 5.94%, MR4 = 5.33%, MR2 = 5.12%).

Sample 5

Code

efa5 <- parameters::factor_analysis(data5, n=5, rotation = "oblimin", sort = TRUE)
plot(efa5)

Code

display(efa5)

Rotated loadings from Factor Analysis (oblimin-rotation)

Variable	MR1	MR4	MR2	MR3	MR5	Complexity	Uniqueness
Burp	0.69	0.06	0.03	-0.10	0.03	1.06	0.50
Wind	0.62	0.01	0.03	-0.02	-0.08	1.04	0.64
Cough	0.54	-0.06	0.01	0.17	0.11	1.31	0.58
Sneeze	0.45	0.02	-0.07	0.21	0.17	1.78	0.61
Bruise	-0.04	0.60	-0.01	0.02	0.06	1.03	0.63
Affective_touch	0.08	0.53	-0.04	-0.08	0.06	1.12	0.67
Blood_Sugar	0.01	0.48	-0.09	0.08	0.19	1.45	0.68
Itch	0.18	0.38	0.08	0.12	-0.02	1.77	0.69
Muscles	0.09	0.38	0.19	3.00e-03	-0.07	1.68	0.77
Thirsty	0.03	-0.17	0.57	-0.01	0.22	1.50	0.55
Temp	2.56e-03	0.18	0.49	0.03	-0.08	1.34	0.71
Hungry	0.02	-0.02	0.48	0.16	0.09	1.32	0.65
Taste	0.19	0.22	0.32	0.05	-0.12	2.84	0.73
Pain	0.10	0.24	0.28	2.96e-03	0.03	2.26	0.77
Breathing	0.04	-0.03	0.03	0.67	-0.07	1.03	0.54
Heart	-0.08	0.05	5.68e-03	0.58	0.06	1.08	0.66
Sex_arousal	0.11	0.14	0.23	0.26	-0.06	3.08	0.76
Urinate	-0.04	0.11	0.11	-0.01	0.61	1.14	0.54
Defecate	0.14	0.03	0.03	-0.02	0.60	1.11	0.53
Vomit	0.15	0.15	-0.04	0.22	0.23	3.56	0.74

The 5 latent factors (oblimin rotation) accounted for 35.25% of the total variance of the original data (MR1 = 8.99%, MR4 = 8.07%, MR2 = 6.47%, MR3 = 5.97%, MR5 = 5.75%).

Sample 6

Code

efa4 <- parameters::factor_analysis(data5, n=4, rotation = "oblimin", sort = TRUE)
plot(efa4)

Code

display(efa4)

Rotated loadings from Factor Analysis (oblimin-rotation)

Variable	MR1	MR4	MR2	MR3	Complexity	Uniqueness
Burp	0.65	0.07	6.06e-03	-0.10	1.07	0.56
Cough	0.63	-0.06	0.03	0.11	1.09	0.56
Sneeze	0.57	0.02	-8.72e-03	0.13	1.10	0.61
Wind	0.54	0.02	-0.03	1.35e-04	1.01	0.71
Defecate	0.34	0.11	0.28	-0.13	2.48	0.70
Vomit	0.27	0.18	0.06	0.13	2.39	0.77
Bruise	-0.06	0.64	-0.04	0.01	1.02	0.63
Affective_touch	0.08	0.55	-0.07	-0.08	1.11	0.67
Blood_Sugar	0.09	0.50	-0.03	0.02	1.07	0.71
Muscles	-2.84e-04	0.41	0.12	0.05	1.19	0.79
Itch	0.15	0.40	0.02	0.13	1.52	0.69
Pain	0.05	0.29	0.26	0.02	2.07	0.78
Taste	0.09	0.24	0.21	0.11	2.76	0.79
Thirsty	0.03	-0.10	0.71	-0.04	1.05	0.51
Hungry	-4.08e-03	0.04	0.51	0.17	1.23	0.65
Temp	-0.09	0.22	0.39	0.10	1.87	0.77
Urinate	0.19	0.18	0.33	-0.11	2.49	0.73
Breathing	0.06	-0.03	-0.01	0.70	1.02	0.50
Heart	-0.01	0.07	0.04	0.50	1.06	0.72
Sex_arousal	0.06	0.17	0.17	0.30	2.35	0.76

The 4 latent factors (oblimin rotation) accounted for 32.02% of the total variance of the original data (MR1 = 9.96%, MR4 = 9.20%, MR2 = 7.36%, MR3 = 5.49%).

Model Construction

Structure

Code

m1 <- "
Interoception =~ Hungry + Thirsty + Urinate + Defecate + Itch + Bruise + Muscles + Pain + Breathing + Heart + Cough + Sneeze + Wind + Burp
"


m7 <- "
# ----
Sustenance =~ Hungry + Thirsty
Elimination =~ Urinate + Defecate
Skin =~ Itch + Bruise
Nociception =~ Muscles + Pain
Affect =~ Breathing + Heart
# ----
Expulsion =~ Cough + Sneeze + Wind + Burp

# --- 
# Ambiguous: Temp + Vomit + Affective_touch + Sex_arousal + Taste
# Discard: Tickle + Blood_Sugar
" 

m8 <- "
# ----
Sustenance =~ Hungry + Thirsty
Elimination =~ Urinate + Defecate
Skin =~ Itch + Bruise
Nociception =~ Muscles + Pain
Affect =~ Breathing + Heart
# ----
Expulsion_Sudden =~ Cough + Sneeze 
Expulsion_Gastric =~ Wind + Burp
" 
# Ambiguous: Temp + Vomit + Affective_touch + Sex_arousal + Taste
# Discard: Tickle + Blood_Sugar

m1_all <- lavaan::cfa(m1, data = data_all)
m7_all <- lavaan::cfa(m7, data = data_all)
m8_all <- lavaan::cfa(m8, data = data_all)

anova(m1_all, m7_all, m8_all) |> 
  parameters::parameters() |> 
  display()

Parameter	AIC	BIC	Chi2	df	p
m8_all	-14218.06	-13919.70	239.58	56
m7_all	-14036.69	-13774.87	432.95	62	< .001
m1_all	-12372.54	-12202.05	2127.10	77	< .001

Anova Table (Type 1 tests)

Higher Order Factors

Code

m8h1 <- paste0(m8, "
              Interoception =~ Sustenance + Elimination + Skin + Nociception + Affect + Expulsion_Sudden + Expulsion_Gastric
              ")

m8h3 <- paste0(m8, "
              Valenced =~ Nociception + Affect + Skin
              Expulsion =~ Expulsion_Sudden + Expulsion_Gastric 
              Homeostasis =~ Sustenance + Elimination
              ")

m8h1_all <- lavaan::cfa(m8h1, data = data_all)
m8h3_all <- lavaan::cfa(m8h3, data = data_all)
anova(m8_all, m8h1_all, m8h3_all) |> 
  parameters::parameters() |> 
  display()

Parameter	AIC	BIC	Chi2	df	p
m8_all	-14218.06	-13919.70	239.58	56
m8h3_all	-14104.40	-13873.02	375.24	67	< .001
m8h1_all	-13903.22	-13690.11	582.42	70	< .001

Anova Table (Type 1 tests)

Model Performance

Code

rbind(
  performance::performance(m8_all, metrics=c("Chi2", "RMSEA", "CFI", "SRMR")) |> 
    mutate(Sample = "All"),
  performance::performance(update(m8_all, data=data1), metrics=c("Chi2", "RMSEA", "CFI", "SRMR")) |> 
    mutate(Sample = "Sample 1"),
  performance::performance(update(m8_all, data=data2), metrics=c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 2"),
  performance::performance(update(m8_all, data=data3), metrics=c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 3"),
  performance::performance(update(m8_all, data=data4), metrics=c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 4"),
  performance::performance(update(m8_all, data=data5), metrics=c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 5"),
  performance::performance(update(m8_all, data=data6), metrics=c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 6"),
  performance::performance(lavaan::cfa(str_remove(m8, fixed("\nExpulsion_Sudden =~ Cough + Sneeze")), data=data7), metrics=c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 7")
) |> 
  display()

Chi2	RMSEA	CFI	SRMR	Sample
239.58	0.03	0.98	0.02	All
88.22	0.04	0.98	0.03	Sample 1
66.50	0.02	0.99	0.03	Sample 2
114.03	0.04	0.98	0.03	Sample 3
115.89	0.03	0.98	0.02	Sample 4
83.61	0.03	0.98	0.03	Sample 5
91.47	0.05	0.97	0.04	Sample 6
82.59	0.05	0.97	0.03	Sample 7

No evidence for higher order factors. The unidimensionality might be a byproduct of the similarities in phrasing of the items.

Study 2 - Correlates

Code

df1 <- cbind(df1, as.data.frame(predict(lavaan::cfa(m8, data = data1)))) |> 
  mutate(Sample = "Sample 1")
df2 <- cbind(df2, as.data.frame(predict(lavaan::cfa(m8, data = data2)))) |> 
  mutate(Sample = "Sample 2")
df3 <- cbind(df3, as.data.frame(predict(lavaan::cfa(m8, data = data3)))) |> 
  mutate(Sample = "Sample 3")
df4 <- cbind(df4, as.data.frame(predict(lavaan::cfa(m8, data = data4)))) |> 
  mutate(Sample = "Sample 4")
df5 <- cbind(df5, as.data.frame(predict(lavaan::cfa(m8, data = data5)))) |> 
  mutate(Sample = "Sample 5")
df6 <- cbind(df6, as.data.frame(predict(lavaan::cfa(m8, data = data6)))) |> 
  mutate(Sample = "Sample 6")
df7 <- cbind(df7, as.data.frame(predict(lavaan::cfa(str_remove(m8, fixed("\nExpulsion_Sudden =~ Cough + Sneeze")), data = data7)))) |> 
  mutate(Sample = "Sample 7", Expulsion_Sudden = NA)

vars <- names(as.data.frame(predict(m8_all)))

Demographics

Code

df_demo <- rbind(
  select(df1, Sample, Age, Gender, all_of(vars)),
  select(df2, Sample, Age, Gender, all_of(vars)),
  select(df3, Sample, Age, Gender, all_of(vars)),
  select(df4, Sample, Age, Gender, all_of(vars)),
  select(df5, Sample, Age, Gender, all_of(vars)),
  select(df6, Sample, Age, Gender, all_of(vars)),
  select(df7, Sample, Age, Gender, all_of(vars))
) 

make_lm_data <- function(df, vars, predictor="Age") {
  dat <- data.frame()
  for(resp in vars) {
      if(length(unique(df$Sample)) == 1) {
        m <- lm(as.formula(paste0(resp, " ~ ", predictor)), data = df)
        param <- parameters::parameters(m)[2, ] |> 
          mutate(Sample = unique(df$Sample)) |> 
          select(-Parameter)
      } else {
        m <- lm(as.formula(paste0(resp, " ~ Sample / ", predictor)), data = df)
        param <- parameters::parameters(m) |> 
          filter(str_detect(Parameter, ":")) |> 
          separate(Parameter, c("Predictor", "Sample"), sep = "Sample ") |> 
          mutate(Sample = paste0("Sample ", str_remove(Sample, fixed(paste0(":", predictor)))))
      }
    
      dat <- param |> 
        mutate(Predictor = predictor,
               Repsonse = resp,
               t_low = CI_low / SE,
               t_high = CI_high / SE) |> 
        rbind(dat)
  }
  dat
}

data <- rbind(
  make_lm_data(df_demo, vars, predictor="Age"),
  make_lm_data(filter(df_demo, Gender %in% c("Male", "Female")), vars, predictor="Gender")
)

BPQ

Code

x <- intersect(names(df3)[str_detect(names(df3), "BPQ")], 
               names(df4)[str_starts(names(df4), "BPQ")])

df_bpq <- rbind(
  select(df3, Sample, all_of(x), all_of(vars)),
  select(df4, Sample, all_of(x), all_of(vars))
) 

# estimate_density(select(df_bpq, Sample, starts_with("BPQ")), at="Sample")  |> 
#   plot() +
#   facet_grid(~Sample)

plot(parameters::n_components(select(df_bpq, starts_with("BPQ"))))

Code

df_bpq$BPQ_PCA1 <- predict(parameters::principal_components(select(df_bpq, starts_with("BPQ")), n=1))$PC1

data <- rbind(data,
              make_lm_data(df_bpq, vars, predictor="BPQ_PCA1"))

# Why opposite effects?

TAS

Code

x <- intersect(names(df3)[str_detect(names(df3), "TAS")], 
               names(df5)[str_starts(names(df5), "TAS")])

PHQ

Code

names(df5)[str_detect(names(df5), "PHQ")]

 [1] "PHQ15_sum" "PHQ9_sum"  "PHQ9_1"    "PHQ9_2"    "PHQ9_3"    "PHQ9_4"   
 [7] "PHQ9_5"    "PHQ9_6"    "PHQ9_7"    "PHQ9_8"    "PHQ15_1"   "PHQ15_2"  
[13] "PHQ15_3"   "PHQ15_4"   "PHQ15_5"   "PHQ15_6"   "PHQ15_7"   "PHQ15_8"  
[19] "PHQ15_9"   "PHQ15_10"  "PHQ15_11"  "PHQ15_12"  "PHQ15_13" 

Code

data <- rbind(data,
              make_lm_data(df5, vars, predictor="PHQ9_sum"),
              make_lm_data(df5, vars, predictor="PHQ15_sum"))

BDI

Code

# names(df1)
# names(df2)
# names(df3)  # TAS, BPQ
# names(df4)  # BPQ
# names(df5)  # PHQ, TAS
# names(df6)  # BDI, STAI
# names(df7)  # MAIA

data <- rbind(data,
              make_lm_data(df6, vars, predictor="BDI2_Total"))

Summary

Code

data |> 
  mutate(sig = ifelse(p < .01, "p < .01", ifelse(p < .001, "p < .001", "N.S."))) |> 
  ggplot(aes(x=Predictor)) +
  geom_hline(yintercept=0, linetype="dashed") +
  geom_pointrange(aes(y=t, ymin=t_low, ymax=t_high, color=sig), position=position_dodge2(width=0.3)) +
  facet_grid(~Repsonse, scales = "free_y") +
  coord_flip()

Discussion

Benefits of the IAS: - Straightforward and sensation-centered items

Recommendations: - Remove Itch (redundant + issue in Chinese) - Use analog scales

Limitations: - Not much clear theorethical or empirical structure (small grouping of items) - Limited variability (clear mode at 4/5) - Ambiguous items which grouping depends on the context (and its awareness) - E.g., heart beating fast, vomit when nauseaous - Few items for some modalities (e.g., 1 for heart) - Positive phrasing: benefits but also might exacerbate positivity bias (and thus unidimensionality)

Need for context-specific items (cross-modal when possible, i.E., cardioception, respiroception, etc.).

New Scale: Multimodal Interoceptive Sensitivity Scale (MISS)